# -*- coding: utf-8 -*-
"""
Created on Mon May 19 23:44:02 2025

@author: 何敏
"""

import pandas as pd
import numpy as np
from linearmodels.panel import PanelOLS, RandomEffects
from scipy import stats

# 假设你的数据框名称是df
df = pd.read_csv("data.csv",encoding='GB2312')

# 设置面板数据格式：省份和年份为索引
df = df.set_index(["省份", "年份"])

# 定义因变量（区域经济发展）和自变量（智能制造水平 + 控制变量）
y = df["区域经济发展"]
X = df[["智能制造水平", "平均受教育年限", "固定资产投资总额", "第三产业增加值比GDP", "一般预算支出比GDP", "城镇化率"]]

# 添加常数项
X = X.assign(const=1)

# 1. 固定效应模型
fixed_effects = PanelOLS(y, X, entity_effects=True)
fixed_effects_results = fixed_effects.fit()

# 2. 随机效应模型
random_effects = RandomEffects(y, X)
random_effects_results = random_effects.fit()

# 3. Hausman 检验
b_fixed = fixed_effects_results.params
b_random = random_effects_results.params

cov_fixed = fixed_effects_results.cov
cov_random = random_effects_results.cov

diff = b_fixed - b_random
stat = np.dot(np.dot(diff.T, np.linalg.inv(cov_fixed - cov_random)), diff)
p_value = stats.chi2.sf(stat, diff.shape[0])  # 计算 p 值

# 输出 Hausman 检验结果
print(f"Hausman 检验统计量: {stat}")
print(f"Hausman 检验 p 值: {p_value}")

if p_value < 0.05:
    print("Hausman 检验显著，选择固定效应模型（FEM）")
else:
    print("Hausman 检验不显著，选择随机效应模型（REM）")